#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gmp.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <pthread.h>
#define BFALSE         0
#define BTRUE          1
#define NUMDIG         100001

/* Define regras para o cálculo do índice cíclico. Por exemplo:
   Isso evita ficar atribuindo de m[1] para m[0] para deixar o valor antigo sempre em m[0] 
   m[0] será o valor antigo quando i for par e, quando ímpar, m[1] será o valor antigo */
#define INDEX(i,j,k) (i+j)%k
#define INDEX2(i,j) INDEX(i,j,2)
#define INDEX3(i,j) INDEX(i,j,3)

mpf_t precisao; //variavel global

/* Método de Gauss-Legendre concorrente */
typedef struct
{
    int i;          //# iterações
    mpf_t a[3];     //variáveis a, b, t, p, pi, para
    mpf_t b[3];
    mpf_t t[2];
    mpf_t p[3];
    mpf_t pi[2];    
    int conv;       /* BTRUE se convergiu */
} parametros; /* variaveis para o algoritmo concorrente */

/* cálculo de a */
void *calcular_a(void *arg)
{
    parametros *pgauss = (parametros *) arg;
    int i0 = INDEX3(pgauss->i, 0);
    int i1 = INDEX3(pgauss->i, 1);

    /* Realiza o cálculo de a */
    mpf_add(pgauss->a[i1], pgauss->a[i0], pgauss->b[i0]); 
    mpf_div_ui(pgauss->a[i1], pgauss->a[i1], 2); // a1 = (a0+b0)/2

    return NULL;
}

/* cálculo de b */
void *calcular_b(void *arg)
{
    parametros *pgauss = (parametros *) arg;
    int i0 = INDEX3(pgauss->i, 0);
    int i1 = INDEX3(pgauss->i, 1);
    
    /* Realiza o cálculo de b */
    mpf_mul(pgauss->b[i1], pgauss->a[i0], pgauss->b[i0]); 
    mpf_sqrt(pgauss->b[i1], pgauss->b[i1]); // b1 = sqrt(a0*b0)

    return NULL;
}

/* cálculo de t */
void *calcular_t(void *arg)
{
    parametros *pgauss = (parametros *) arg;
    int i0, i1, j0, j1;

    /* "Bolhas" para cálculo de t */
    if(pgauss->i < 1) return NULL;

    i0 = INDEX3(pgauss->i-1, 0);
    i1 = INDEX3(pgauss->i-1, 1);
    j0 = INDEX2(pgauss->i-1, 0);
    j1 = INDEX2(pgauss->i-1, 1);

    /* cálculo de t */
    mpf_sub(pgauss->t[j1], pgauss->a[i0], pgauss->a[i1]); 
    mpf_mul(pgauss->t[j1], pgauss->t[j1], pgauss->t[j1]); 
    mpf_mul(pgauss->t[j1], pgauss->t[j1], pgauss->p[i0]); 
    mpf_sub(pgauss->t[j1], pgauss->t[j0], pgauss->t[j1]); // t1 = t0-p0*(a0-a1)^2

    return NULL;
}

/* cálculo de p */
void *calcular_p(void *arg)
{
    parametros *pgauss = (parametros *) arg;
    int i0 = INDEX3(pgauss->i, 0);
    int i1 = INDEX3(pgauss->i, 1);

    /* cálculo de p */
    mpf_mul_ui(pgauss->p[i1], pgauss->p[i0], 2); // p1 = 2*p0

    return NULL;
}

/* cálculo de Pi */
void *calcular_pi(void *arg)
{
    parametros *pgauss = (parametros *) arg;
    int j0, i1, j1;

    /* "Bolhas" para cálculo do Pi */
    if(pgauss->i < 2) return NULL;

    i1 = INDEX3(pgauss->i-2, 1);
    j0 = INDEX2(pgauss->i-2, 0);
    j1 = INDEX2(pgauss->i-2, 1);

    /* cálculo de Pi */
    mpf_add(pgauss->pi[j1], pgauss->a[i1], pgauss->b[i1]); 
    mpf_mul(pgauss->pi[j1], pgauss->pi[j1], pgauss->pi[j1]); 
    mpf_div(pgauss->pi[j1], pgauss->pi[j1], pgauss->t[j1]); 
    mpf_div_ui(pgauss->pi[j1], pgauss->pi[j1], 4); // pi = (a1+b1)^2/(4*t1)

    /* verifica convergencia */
    mpf_t aux;
    int convergencia;
    mpf_init(aux);

    mpf_sub(aux, pgauss->pi[j1], pgauss->pi[j0]); 
    mpf_abs(aux, aux);                            // aux = |pgauss->pi[j1] - pgauss->pi[j0]|
    convergencia = (mpf_cmp(aux, precisao) <= 0);
    mpf_clear(aux);

    if(convergencia) 
      pgauss->conv = BTRUE;

    return NULL;
}

int main()
{
    mpf_t pi;
    int i,j;
    //VARIAVEIS PARA GUARDAR O TEMPO
    int tempo;
    struct timeval inicio, fim;
    
    gmp_printf("calculando precisao... \n");
    /* Precisão em bits = ceil( 10000001 * log2(10) ) */
    //mpf_set_default_prec(33219281);
    mpf_set_default_prec(33219281);
    mpf_init_set_ui(precisao, 10);
    mpf_pow_ui(precisao, precisao, NUMDIG); 
    mpf_ui_div(precisao, 1, precisao); // precisao = 1/(10^PREC)
    //gmp_printf("precisao = %F\n\n", precisao);
    mpf_init(pi);    

    pthread_t thread[5];
    parametros pgauss;
    
    /* Inicializa estrutura paralela */

    /* Valores iniciais para o método de Gauss-Legendre */
    /* PEGA O TEMPO INICIAL DE EXECUCAO */
    gettimeofday(&inicio, NULL);
    
    mpf_init_set_ui(pgauss.a[0], 1);           // a0 = 1
    mpf_init_set_ui(pgauss.b[0], 2);
    mpf_sqrt(pgauss.b[0], pgauss.b[0]);
    mpf_ui_div(pgauss.b[0], 1, pgauss.b[0]);   // b0 = 1/sqrt(2)
    mpf_init_set_d(pgauss.t[0], 0.25);         // t0 = 1/4
    mpf_init_set_ui(pgauss.p[0], 1);           // p0 = 1

    mpf_inits(pgauss.a[1], pgauss.a[2], pgauss.b[1], pgauss.b[2], pgauss.t[1], pgauss.p[1], pgauss.p[2], pgauss.pi[0], pgauss.pi[1], NULL);

    pgauss.conv = BFALSE;
    for(i = 0; !pgauss.conv; i++)
    {
      pgauss.i = i;

      /* Cria uma thread para cada parâmetro */
      pthread_create(&thread[0], NULL, calcular_a, &pgauss);
      pthread_create(&thread[1], NULL, calcular_b, &pgauss);
      pthread_create(&thread[2], NULL, calcular_t, &pgauss);
      pthread_create(&thread[3], NULL, calcular_p, &pgauss);
      pthread_create(&thread[4], NULL, calcular_pi, &pgauss);
      for(j = 0; j < 5; j++) pthread_join(thread[j], NULL);
    }

    /* pi calculado */
    mpf_set(pi, pgauss.pi[INDEX2(i, 1)]);
    gmp_printf("numero de iteracoes: %i \nvalor estimado de pi: %.1000Ff \n",i,pi);
    
     //PEGA O TEMPO FINAL DE EXECUCAO
    gettimeofday(&fim, NULL);

    //DIFERENCA ENTRE FINAL E INICIAL
    tempo = fim.tv_sec - inicio.tv_sec;
    printf("TEMPO DE EXECUÇÃO %d segundos\n", tempo);

    gmp_printf("imprimindo arquivo 'pi-gauss-legendre-con.txt'\n");
    FILE *fp = fopen("pi-gauss-legendre-con.txt", "w"); /* Abre o arquivo */
    mpf_out_str(fp, 10, NUMDIG, pi); /* Imprime com precisão 10^7 dígitos*/   
    gmp_printf("\n");
    fclose(fp); 
    return 0;
    
    /* clear nos parametros */
    mpf_clears(pgauss.a[0], pgauss.a[1], pgauss.a[2], pgauss.b[0], pgauss.b[1], pgauss.b[2], pgauss.t[0], pgauss.t[1],
            pgauss.p[0], pgauss.p[1], pgauss.p[2], pgauss.pi[0], pgauss.pi[1], NULL);

    mpf_clear(precisao);
    mpf_clear(pi);
    return 0;
}

